!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

 module diagnostics

!BOP
! !MODULE: diagnostics
! !DESCRIPTION:
!  This module contains routines and data for producing diagnostic
!  output, including global diagnostics, cfl checks and transport
!  diagnostics.
!
! !REVISION HISTORY:
!  SVN:$Id: diagnostics.F90 28439 2011-05-18 21:40:58Z njn01 $
!
! !USES:

   use POP_KindsMod
   use POP_IOUnitsMod
   use POP_SolversMod

   use domain
   use constants
   use prognostic
   use time_management
   use io
   use broadcast
   use global_reductions
   use grid
   use forcing
   use forcing_fields
   use timers
   use exit_mod
   use tavg, only: define_tavg_field,     &
                   accumulate_tavg_field, accumulate_tavg_now, &
                   tavg_method_avg, tavg_method_max, tavg_method_min
   use movie, only: define_movie_field, movie_requested, update_movie_field
   use vmix_kpp, only: HMXL, KPP_HBLT
   use registry
   use io_tools
   use gather_scatter

   implicit none
   private
   save

! !PUBLIC MEMBER FUNCTIONS:

   public :: init_diagnostics,        &
             diag_init_sums,          &
             diag_global_preupdate,   &
             diag_global_afterupdate, &
             diag_print,              &
             diag_transport,          &
             cfl_advect,              &
             cfl_vdiff,               &
             cfl_hdiff,               &
             cfl_check,               &
             check_KE,                &
             diag_velocity

! !PUBLIC DATA MEMBERS:

   logical (log_kind), public :: &
      ldiag_global,          &! time to compute global diagnostics
      ldiag_cfl,             &! time to compute cfl diagnostics
      ldiag_transport,       &! time to compute transport diagnostics
      ldiag_velocity,        &! compute velocity diagnostics
      cfl_all_levels,        &! writes cfl  diags for all vert levels
      diag_all_levels         ! writes some diags for all vert levels

   !***  arrays for holding various diagnostic results
   !***  public for now as they are modified directly in baroclinic

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
      public :: &
      DIAG_KE_ADV_2D,   &!
      DIAG_KE_HMIX_2D,  &!
      DIAG_KE_VMIX_2D,  &!
      DIAG_KE_PRESS_2D, &!
      DIAG_PE_2D

   real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), &
      public :: &
      DIAG_TRACER_ADV_2D,    &!
      DIAG_TRACER_HDIFF_2D,  &!
      DIAG_TRACER_VDIFF_2D,  &!
      DIAG_TRACER_SOURCE_2D, &!
      DIAG_TRACER_SFC_FLX

!-----------------------------------------------------------------------
!
!     global budget diagnostics
!
!-----------------------------------------------------------------------

   real (r8), dimension(nt),public ::   &
     tracer_mean_initial,               &! SUM [volume*tracer] at the beginning 
                                         !    of a time-averaging interval
     tracer_mean_final                   ! SUM [volume*tracer] at the end of a
                                         ! time-averaging interval

   real (r8), public                ::  &
      volume_t_initial,                 & ! T-point volume at the beginning of a
                                          ! time-averaging interval
      volume_t_final                      ! T-point volume at the end of a
                                          ! time-averaging interval

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  variables for controlling frequency and output of diagnostics
!
!-----------------------------------------------------------------------

   integer (int_kind) ::     &
      diag_unit,             &! i/o unit for output diagnostic file
      trans_unit,            &! i/o unit for output transport file
      velocity_unit           ! i/o unit for output velocity file

   integer (int_kind) ::     &
      diag_global_flag,      &! time flag id for global diags
      diag_cfl_flag,         &! time flag id for cfl    diags
      diag_transp_flag        ! time flag id for transp diags

   character (char_len) ::   &
      diag_outfile,                  &! current  filename for diagnostic output
      diag_transport_outfile,        &! current  filename for transport output
      diag_velocity_outfile           ! current  filename for velocity output

!-----------------------------------------------------------------------
!
!  variables specific to ccsm coupled code
!
!-----------------------------------------------------------------------
   logical (log_kind) :: &
      lccsm = .false.

   character (char_len)     ::  &
      ccsm_diag_date

!-----------------------------------------------------------------------
!
!  global diagnostics
!
!-----------------------------------------------------------------------

   real (r8) ::        &
      diag_ke,         &! mean KE at new time
      diag_ke_adv,     &! KE change due to advection
      diag_ke_free,    &! advective KE change due to free surface
      diag_ke_hmix,    &! KE change due to horizontal diffusion
      diag_ke_vmix,    &! KE change due to vertical   diffusion
      diag_ke_press,   &! KE change due to pressure gradient
      diag_pe,         &! change in potential energy
      diag_press_free, &! pressure work contribution from free surface
      diag_ke_psfc,    &! KE change due to surface pressure
      diag_ws,         &! KE change due to wind stress
      diag_sealevel     ! global mean seal level

   real (r8), dimension(nt) :: &
      dtracer_avg,             &! avg tracer change for this time step
      dtracer_abs,             &! abs tracer change for this time step
      diag_tracer_adv,         &! tracer change due to advection
      diag_tracer_hdiff,       &! tracer change due to horiz diffusion
      diag_tracer_vdiff,       &! tracer change due to vert  diffusion
      diag_tracer_source,      &! tracer change due to source terms
      avg_tracer,              &! global average tracer at new time
      sfc_tracer_flux           ! sfc tracer flux (Fw*tracer) or
                                !   resid sfc flux (w*tracer)

   real (r8), dimension(km,nt) :: &
      avg_tracer_k                ! mean tracer at every level


!-----------------------------------------------------------------------
!
!  cfl diagnostics
!
!-----------------------------------------------------------------------

   real (r8), dimension(max_blocks_clinic,km) :: &
      cfl_advuk_block,     &! zonal advective cfl number for level k
      cfl_advvk_block,     &! merid advective cfl number for level k
      cfl_advwk_block,     &! vert  advective cfl number for level k
      cfl_advtk_block,     &! total advective cfl number for level k
      cfl_vdifftk_block,   &! tracer   vert  diff cfl num for level k
      cfl_vdiffuk_block,   &! momentum vert  diff cfl num for level k
      cfl_hdifftk_block,   &! tracer   horiz diff cfl num for level k
      cfl_hdiffuk_block     ! momentum horiz diff cfl num for level k

   integer (int_kind), dimension(2,max_blocks_clinic,km) :: &
      cfladd_advuk_block,       &! horiz addresses for max cfl numbers
      cfladd_advvk_block,       &!   at level k
      cfladd_advwk_block,       &
      cfladd_advtk_block,       &
      cfladd_vdifftk_block,     &
      cfladd_vdiffuk_block,     &
      cfladd_hdifftk_block,     &
      cfladd_hdiffuk_block

!-----------------------------------------------------------------------
!
!  transport diagnostics
!
!-----------------------------------------------------------------------

   type :: transport
      character(char_len) :: name   ! name for transport
      integer (int_kind)  :: type   ! type (meridional or zonal)
      integer (int_kind), dimension(6,max_blocks_clinic) :: &
                             add    ! address range for computing transp
      real (r8)           :: mass, &! mass transport
                             heat, &! heat transport
                             salt   ! salt transport
   end type

   integer (int_kind), parameter :: &
      transport_type_zonal = 1, &! zonal or meridional transport type
      transport_type_merid = 2

   integer (int_kind) :: &
      num_transports      ! number of transport diagnostics computed

   type (transport), dimension(:), allocatable :: &
      transports          ! transport info for all requested transports

!-----------------------------------------------------------------------
!
!  ids for tavg diagnostics computed from diagnostic
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_HMXL,         &! tavg id for average mixed layer depth
      tavg_HMXL_2,       &! tavg id for average mixed layer depth, stream #2 (allows two frequencies)
      tavg_XMXL,         &! tavg id for maximum mixed layer depth
      tavg_XMXL_2,       &! tavg id for maximum mixed layer depth, stream #2
      tavg_TMXL,         &! tavg id for minimum mixed layer depth
      tavg_HBLT,         &! tavg id for average boundary layer depth
      tavg_XBLT,         &! tavg id for maximum boundary layer depth
      tavg_TBLT           ! tavg id for minimum boundary layer depth

!-----------------------------------------------------------------------
!
!  movie ids 
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      movie_HMXL,          &! movie id for average mixed layer depth
      movie_HBLT            ! movie id for average boundary layer depth

!-----------------------------------------------------------------------
!
!  solver diagnostics
!
!-----------------------------------------------------------------------

   integer (POP_i4) :: &
      iterationCount    ! current solver iteration count

   real (POP_r8) ::    &
      rmsResidual       ! mean solver residual

!-----------------------------------------------------------------------
!
!  velocity diagnostics
!
!-----------------------------------------------------------------------
 
   integer (int_kind), parameter ::  &
      num_vel_loc = 5
 
   integer (int_kind), dimension (num_vel_loc) ::     &
      i_loc, j_loc
 
   real (r8), dimension (num_vel_loc) ::     &
      vlat_loc, vlon_loc 
 


!EOC
!***********************************************************************

 contains

!***********************************************************************
!BOP
! !IROUTINE: init_diagnostics
! !INTERFACE:

 subroutine init_diagnostics

! !DESCRIPTION:
!  Initializes diagnostic module quantities.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      nu,                &! unit for contents input file
      iblock, n,         &! dummy loop indices
      nloc, nreach,      &!   ditto
      nml_error,         &! namelist i/o error flag
      ib,ie,jb,je         ! beg,end indices for block physical domain

   integer (int_kind), dimension(6) :: &
      tmp_add             ! transport global addresses
                          ! (ibeg,iend,jbeg,jend,kbeg,kend)

   character (char_len) ::  &
      diag_global_freq_opt, &! choice for freq of global diagnostics
      diag_cfl_freq_opt,    &! choice for freq of cfl diagnostics
      diag_transp_freq_opt, &! choice for freq of transport diagnostics
      diag_transport_file,  &! filename for choosing fields for output
      transport_ctype,      &! type of transport (zonal,merid)
      outfile_tmp,          &! temp for appending to outfile name
      string                 ! temp for creating outfile name

   integer (int_kind) ::     &
      diag_global_freq_iopt, &! freq option for computing global diags
      diag_global_freq,      &! freq for computing global diagnostics
      diag_cfl_freq_iopt,    &! freq option for computing cfl diags
      diag_cfl_freq,         &! freq for computing cfl diagnostics
      diag_transp_freq_iopt, &! freq option for comp transport diags
      diag_transp_freq        ! freq for computing transp diagnostics

   character (10) ::        &
      cdate                 ! character date

   character (47), parameter :: &! output formats for freq options
      gfreq_fmt = "('Global    diagnostics computed every ',i8,a8)", &
      cfreq_fmt = "('CFL       diagnostics computed every ',i8,a8)", &
      tfreq_fmt = "('Transport diagnostics computed every ',i8,a8)"

   namelist /diagnostics_nml/diag_global_freq_opt, diag_global_freq, &
                             diag_cfl_freq_opt,    diag_cfl_freq,    &
                             diag_transp_freq_opt, diag_transp_freq, &
                             diag_transport_file,                    &
                             diag_all_levels, cfl_all_levels,        &
                             diag_outfile, diag_transport_outfile,   &
                             diag_velocity_outfile,                  &
                             ldiag_velocity

   type (block) ::         &
      this_block           ! block information for current block

!-----------------------------------------------------------------------
!
!  local variables for velocity diagnostics lat/lon search
!
!-----------------------------------------------------------------------

   integer (int_kind), parameter :: search_error = -999

   real (r8)     :: &
      diag_vel_search_reach,  &! size of lat/lon search length
      diag_vel_search_inc,    &! increment in " "
      lat_reach,              &! local search-length names
      lon_reach,              &
      min_distance,           &!smallest distance from desired lat/lon point
      distance                 !actual distance from desired lat/lon point

   integer (int_kind) ::      &
       i,j,                   &! dummy indices
       diag_vel_search_max

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: WORK

   real (r8), allocatable, dimension (:,:) :: &
      ULAT_G,ULON_G   

   integer (int_kind), allocatable, dimension(:,:) ::  &
      REGION_MASK_G

   real (r8), parameter, dimension (num_vel_loc) ::     &
      lon_loc = (/ 220.0_r8, 220.0_r8, 335.0_r8, 335.0_r8, 250.0_r8 /)
 
   real (r8), parameter, dimension (num_vel_loc) ::     &
      lat_loc = (/   0.0_r8,   2.0_r8,   0.0_r8,   2.0_r8,   0.0_r8 /)

!-----------------------------------------------------------------------
!
!  determine if this is a ccsm coupled run
!
!-----------------------------------------------------------------------

   lccsm = registry_match('lccsm')


!-----------------------------------------------------------------------
!
!  read diagnostic file output frequency and filenames from namelist
!
!-----------------------------------------------------------------------

   !*** set default values

   diag_global_freq_opt   = 'never'
   diag_cfl_freq_opt      = 'never'
   diag_transp_freq_opt   = 'never'
   diag_transport_file    = 'unknown_transport_file'
   diag_outfile           = 'unknown_diag_outfile'
   diag_transport_outfile = 'unknown_transport_outfile'
   diag_all_levels        = .false.
   cfl_all_levels         = .false.
   ldiag_velocity         = .false.
   diag_velocity_outfile  = 'unknown_velocity_outfile'

   if (my_task == master_task) then
      open (nml_in, file=nml_filename, status='old',iostat=nml_error)
      if (nml_error /= 0) then
         nml_error = -1
      else
         nml_error =  1
      endif
      do while (nml_error > 0)
         read(nml_in, nml=diagnostics_nml,iostat=nml_error)
      end do
      if (nml_error == 0) close(nml_in)
   endif

   call broadcast_scalar(nml_error, master_task)
   if (nml_error /= 0) then
      call exit_POP(sigAbort,'ERROR reading diagnostics_nml')
   endif

!-----------------------------------------------------------------------
!  ccsm filenames must meet ccsm output-file naming conventions
!  actual ccsm filenames will be constructed in diag_print
!-----------------------------------------------------------------------
   if (lccsm) then

     call ccsm_date_stamp (ccsm_diag_date, 'ymds')

     string = diag_outfile
     diag_outfile = trim(string)/&
                                 &/'.'/&
                                 &/trim(ccsm_diag_date)

     string = diag_transport_outfile
     diag_transport_outfile = trim(string)/&
                                           &/'.'/&
                                           &/trim(ccsm_diag_date)

     string = diag_velocity_outfile
     diag_velocity_outfile  = trim(string)/&
                                           &/'.'/&
                                           &/trim(ccsm_diag_date)

   else

!-----------------------------------------------------------------------
!  non-ccsm filenames:
!  append runid, initial date to output file names
!  concatenation operator must be split across lines to avoid problems
!    with preprocessors
!
!-----------------------------------------------------------------------

     if (date_separator == ' ') then
        cdate(1:4) = cyear
        cdate(5:6) = cmonth
        cdate(7:8) = cday
        cdate(9:10)= '  '
     else
        cdate(1:4) = cyear
        cdate(5:5) = date_separator
        cdate(6:7) = cmonth
        cdate(8:8) = date_separator
        cdate(9:10) = cday
     endif

     outfile_tmp = char_blank
     outfile_tmp = trim(diag_outfile)/&
                                      &/'.'/&
                                      &/trim(runid)/&
                                      &/'.'/&
                                      &/trim(cdate)
     diag_outfile = char_blank
     diag_outfile = trim(outfile_tmp)
  
     outfile_tmp = char_blank
     outfile_tmp = trim(diag_transport_outfile)/&
                                                &/'.'/&
                                                &/trim(runid)/&
                                                &/'.'/&
                                                &/trim(cdate)
     diag_transport_outfile = char_blank
     diag_transport_outfile = trim(outfile_tmp)

   endif ! lccsm

!-----------------------------------------------------------------------
!
!  set, broadcast and print options to log file
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(stdout,blank_fmt)
      write(stdout,delim_fmt)
      write(stdout,blank_fmt)
      write(stdout,'(a18)') 'Diagnostic options'
      write(stdout,blank_fmt)
      write(stdout,delim_fmt)

      select case (diag_global_freq_opt)
      case ('never')
         diag_global_freq_iopt = freq_opt_never
         write(stdout,'(a31)') 'Global diagnostics not computed'
      case ('nyear')
         diag_global_freq_iopt = freq_opt_nyear
         write(stdout,gfreq_fmt) diag_global_freq, ' years  '
      case ('nmonth')
         diag_global_freq_iopt = freq_opt_nmonth
         write(stdout,gfreq_fmt) diag_global_freq, ' months '
      case ('nday')
         diag_global_freq_iopt = freq_opt_nday
         write(stdout,gfreq_fmt) diag_global_freq, ' days   '
      case ('nhour')
         diag_global_freq_iopt = freq_opt_nhour
         write(stdout,gfreq_fmt) diag_global_freq, ' hours  '
      case ('nsecond')
         diag_global_freq_iopt = freq_opt_nsecond
         write(stdout,gfreq_fmt) diag_global_freq, ' seconds'
      case ('nstep')
         diag_global_freq_iopt = freq_opt_nstep
         write(stdout,gfreq_fmt) diag_global_freq, ' steps  '
      case default
         diag_global_freq_iopt = -1000
      end select

      if (diag_global_freq_iopt /= freq_opt_never) then
           write(stdout,'(a36,a)') &
              'Global diagnostics written to file: ', trim(diag_outfile)
         if (diag_all_levels)     &
            write(stdout,'(a42)') &
               'Diagnostics output for all vertical levels'
      endif

      if (lccsm) then
        write (stdout,'(a39,a)') 'Equatorial velocities written to file: ', &
            trim(diag_velocity_outfile)
      endif

      select case (diag_cfl_freq_opt)
      case ('never')
         diag_cfl_freq_iopt = freq_opt_never
         write(stdout,'(a28)') 'CFL diagnostics not computed'
      case ('nyear')
         diag_cfl_freq_iopt = freq_opt_nyear
         write(stdout,cfreq_fmt) diag_cfl_freq, ' years  '
      case ('nmonth')
         diag_cfl_freq_iopt = freq_opt_nmonth
         write(stdout,cfreq_fmt) diag_cfl_freq, ' months '
      case ('nday')
         diag_cfl_freq_iopt = freq_opt_nday
         write(stdout,cfreq_fmt) diag_cfl_freq, ' days   '
      case ('nhour')
         diag_cfl_freq_iopt = freq_opt_nhour
         write(stdout,cfreq_fmt) diag_cfl_freq, ' hours  '
      case ('nsecond')
         diag_cfl_freq_iopt = freq_opt_nsecond
         write(stdout,cfreq_fmt) diag_cfl_freq, ' seconds'
      case ('nstep')
         diag_cfl_freq_iopt = freq_opt_nstep
         write(stdout,cfreq_fmt) diag_cfl_freq, ' steps  '
      case default
         diag_cfl_freq_iopt = -1000
      end select

      if (diag_cfl_freq_iopt /= freq_opt_never) then
           write(stdout,'(a33,a)') &
              'CFL diagnostics written to file: ',trim(diag_outfile)
         if (cfl_all_levels) then
            write(stdout,'(a46)') &
               'CFL diagnostics output for all vertical levels'
         endif
      endif

      select case (diag_transp_freq_opt)
      case ('never')
         diag_transp_freq_iopt = freq_opt_never
         write(stdout,'(a34)') 'Transport diagnostics not computed'
      case ('nyear')
         diag_transp_freq_iopt = freq_opt_nyear
         write(stdout,tfreq_fmt) diag_transp_freq, ' years  '
      case ('nmonth')
         diag_transp_freq_iopt = freq_opt_nmonth
         write(stdout,tfreq_fmt) diag_transp_freq, ' months '
      case ('nday')
         diag_transp_freq_iopt = freq_opt_nday
         write(stdout,tfreq_fmt) diag_transp_freq, ' days   '
      case ('nhour')
         diag_transp_freq_iopt = freq_opt_nhour
         write(stdout,tfreq_fmt) diag_transp_freq, ' hours  '
      case ('nsecond')
         diag_transp_freq_iopt = freq_opt_nsecond
         write(stdout,tfreq_fmt) diag_transp_freq, ' seconds'
      case ('nstep')
         diag_transp_freq_iopt = freq_opt_nstep
         write(stdout,tfreq_fmt) diag_transp_freq, ' steps  '
      case default
         diag_transp_freq_iopt = -1000
      end select

   endif

   call broadcast_scalar(diag_global_freq_iopt, master_task)
   call broadcast_scalar(diag_global_freq,      master_task)
   call broadcast_scalar(diag_cfl_freq_iopt,    master_task)
   call broadcast_scalar(diag_cfl_freq,         master_task)
   call broadcast_scalar(diag_transp_freq_iopt, master_task)
   call broadcast_scalar(diag_transp_freq,      master_task)
   call broadcast_scalar(diag_all_levels,       master_task)
   call broadcast_scalar(cfl_all_levels,        master_task)
   call broadcast_scalar(ldiag_velocity,        master_task)

   if (diag_global_freq_iopt == -1000) then
      call exit_POP(sigAbort, &
                    'ERROR: unknown global diag frequency option')
   endif
   if (diag_cfl_freq_iopt == -1000) then
      call exit_POP(sigAbort, &
                    'ERROR: unknown cfl diagnostic frequency option')
   endif
   if (diag_transp_freq_iopt == -1000) then
      call exit_POP(sigAbort, &
                    'ERROR: unknown transport diag frequency option')
   endif


!-----------------------------------------------------------------------
!
!  create time flags for computing diagnostics
!
!-----------------------------------------------------------------------

   call init_time_flag('diag_global', diag_global_flag, default=.false.,  &
                                      owner    = 'init_diagnostics',      &
                                      freq_opt = diag_global_freq_iopt,   &
                                      freq     = diag_global_freq)

   call init_time_flag('diag_cfl',    diag_cfl_flag, default=.false.,     &
                                      owner    = 'init_diagnostics',      &
                                      freq_opt = diag_cfl_freq_iopt,      &
                                      freq     = diag_cfl_freq)

   call init_time_flag('diag_transp', diag_transp_flag, default=.false.,  &
                                      owner    = 'init_diagnostics',      &
                                      freq_opt = diag_transp_freq_iopt,   &
                                      freq     = diag_transp_freq)

!-----------------------------------------------------------------------
!
!  read transports file to determine which transports to compute
!
!-----------------------------------------------------------------------

   if (diag_transp_freq_iopt /= freq_opt_never) then

      call get_unit(nu)
      if (my_task == master_task) then
         write(stdout,'(a39,a)') &
            'Transport diagnostics written to file: ', &
             trim(diag_transport_outfile)
         write(stdout,blank_fmt)

         open(nu, file=diag_transport_file, status='old')
         read(nu,*) num_transports

         write(stdout,'(a14,i4,1x,a27)') 'The following ',num_transports, &
                                    'transports will be computed'
      endif

      call broadcast_scalar(num_transports, master_task)

      allocate( transports(num_transports) )

      do n=1,num_transports
         if (my_task == master_task) then
            read(nu,'(6(i4,1x),a5,2x,a)') tmp_add, &
                                 transport_ctype, transports(n)%name
            write(stdout,'(a2,a)') '  ',trim(transports(n)%name)

            select case (transport_ctype(1:5))
            case ('zonal')
               transports(n)%type = transport_type_zonal
            case ('merid')
               transports(n)%type = transport_type_merid
            case default
               transports(n)%type = -1000
            end select
         endif

         call broadcast_scalar(transports(n)%type,master_task)
         call broadcast_array(tmp_add,master_task)
         if (transports(n)%type == -1000) then
            call exit_POP(sigAbort,'ERROR: unknown transport type')
         endif

         !***
         !*** compute local addresses for transport
         !***

         !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je)

         do iblock=1,nblocks_clinic

            this_block = get_block(blocks_clinic(iblock),iblock)  

            ib = this_block%ib
            ie = this_block%ie
            jb = this_block%jb
            je = this_block%je

            transports(n)%add(1,iblock) = nx_block
            transports(n)%add(2,iblock) = 1
            transports(n)%add(3,iblock) = ny_block
            transports(n)%add(4,iblock) = 1
            transports(n)%add(5,iblock) = tmp_add(5)  ! beg k index
            transports(n)%add(6,iblock) = tmp_add(6)  ! end k index

            if (tmp_add(1) <= this_block%i_glob(ie) .and. &
                tmp_add(2) >= this_block%i_glob(ib) ) then

               if (tmp_add(1) >= this_block%i_glob(ib)) then
                  transports(n)%add(1,iblock) = &
                              tmp_add(1) - this_block%i_glob(ib) + ib
               else
                  transports(n)%add(1,iblock) = ib
               endif
               if (tmp_add(2) <= this_block%i_glob(ie)) then
                  transports(n)%add(2,iblock) = &
                              tmp_add(2) - this_block%i_glob(ib) + ib
               else
                  transports(n)%add(2,iblock) = ie
               endif
            endif

            if (tmp_add(3) <= this_block%j_glob(je) .and. &
                tmp_add(4) >= this_block%j_glob(jb) ) then
               if (tmp_add(3) >= this_block%j_glob(jb)) then
                  transports(n)%add(3,iblock) = &
                              tmp_add(3) - this_block%j_glob(jb) + jb
               else
                  transports(n)%add(3,iblock) = jb
               endif
              if (tmp_add(4) <= this_block%j_glob(je)) then
                 transports(n)%add(4,iblock) = &
                              tmp_add(4) - this_block%j_glob(jb) + jb
              else
                 transports(n)%add(4,iblock) = je
              endif
            endif

         end do ! block loop

         !$OMP END PARALLEL DO

      end do

      if (my_task == master_task) close(nu)
      call release_unit(nu)
   endif

!-----------------------------------------------------------------------
!
!  open output files if required, then write a blank and close.
!
!-----------------------------------------------------------------------

   if (diag_global_freq_iopt /= freq_opt_never .or. &
       diag_cfl_freq_iopt    /= freq_opt_never) then

      call get_unit(diag_unit)
      if (my_task == master_task) then
          open(diag_unit, file=diag_outfile, status='unknown')
          write(diag_unit,*)' '
          close(diag_unit)
      endif
   endif

   if (diag_transp_freq_iopt /= freq_opt_never) then
      call get_unit(trans_unit)
      if (my_task == master_task) then
         open(trans_unit, file=diag_transport_outfile,status='unknown')
         write(trans_unit,*)' '
         close(trans_unit)
      endif
   endif

   if ( lccsm ) then
     call get_unit(velocity_unit)
     if (my_task == master_task) then
       open(velocity_unit, file=diag_velocity_outfile, status='unknown')
       write(velocity_unit,*)' '
       close(velocity_unit)
     endif
   endif

!-----------------------------------------------------------------------
!
!  initialize cfl arrays if required
!
!-----------------------------------------------------------------------

   if (diag_cfl_freq_iopt /= freq_opt_never) then
      cfl_advuk_block       = c0
      cfl_advvk_block       = c0
      cfl_advwk_block       = c0
      cfl_advtk_block       = c0
      cfl_vdifftk_block     = c0
      cfl_vdiffuk_block     = c0
      cfl_hdifftk_block     = c0
      cfl_hdiffuk_block     = c0
      cfladd_advuk_block    = 0
      cfladd_advvk_block    = 0
      cfladd_advwk_block    = 0
      cfladd_advtk_block    = 0
      cfladd_vdifftk_block  = 0
      cfladd_vdiffuk_block  = 0
      cfladd_hdifftk_block  = 0
      cfladd_hdiffuk_block  = 0
   endif

!-----------------------------------------------------------------------
!
!  define tavg diagnostic fields
!
!-----------------------------------------------------------------------

   call define_tavg_field(tavg_HMXL,'HMXL',2,                         &
                          tavg_method=tavg_method_avg,                &
                          long_name='Mixed-Layer Depth',              &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_HMXL_2,'HMXL_2',2,                     &
                          tavg_method=tavg_method_avg,                &
                          long_name='Mixed-Layer Depth',              &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_XMXL,'XMXL',2,                         &
                          tavg_method=tavg_method_max,                &
                          long_name='Maximum Mixed-Layer Depth',      &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_XMXL_2,'XMXL_2',2,                     &
                          tavg_method=tavg_method_max,                &
                          long_name='Maximum Mixed-Layer Depth',      &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_TMXL,'TMXL',2,                         &
                          tavg_method=tavg_method_min,                &
                          long_name='Minimum Mixed-Layer Depth',      &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_HBLT,'HBLT',2,                         &
                          tavg_method=tavg_method_avg,                &
                          long_name='Boundary-Layer Depth',           &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_XBLT,'XBLT',2,                     &
                          tavg_method=tavg_method_max,                &
                          long_name='Maximum Boundary-Layer Depth',   &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_TBLT,'TBLT',2,                     &
                          tavg_method=tavg_method_min,                &
                          long_name='Minimum Boundary-Layer Depth',   &
                          units='centimeter', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

!-----------------------------------------------------------------------
!
!  define movie diagnostic fields
!
!-----------------------------------------------------------------------

   call define_movie_field(movie_HMXL,'HMXL',0,                       &
                          long_name='Mixed-Layer Depth',              &
                          units='cm', grid_loc='2110')

   call define_movie_field(movie_HBLT,'HBLT',0,                       &
                          long_name='Boundary-Layer Depth',           &
                          units='cm', grid_loc='2110')

!-----------------------------------------------------------------------
!
!  find global i,j indices close to specified lat,lon points
!
!-----------------------------------------------------------------------
 
  if (ldiag_velocity) then

   j_loc = search_error
   i_loc = search_error
 
   allocate (ULAT_G(nx_global,ny_global), ULON_G(nx_global,ny_global))
   allocate (REGION_MASK_G(nx_global,ny_global))
 
   WORK=ULON*radian
   call gather_global(ULON_G,WORK,master_task,distrb_clinic)
  
   WORK=ULAT*radian
   call gather_global(ULAT_G,WORK,master_task,distrb_clinic)
 
   call gather_global(REGION_MASK_G,REGION_MASK,master_task,distrb_clinic)
 
   diag_vel_search_reach = c0
   diag_vel_search_inc   = 0.1_r8
   diag_vel_search_max   = 10
  
   if (my_task == master_task) then
 
   do nloc = 1, num_vel_loc

     min_distance  = c10
     lat_reach = diag_vel_search_reach 
     lon_reach = diag_vel_search_reach 
 
     reach_loop: do nreach = 1,diag_vel_search_max
 
       lat_reach = lat_reach + diag_vel_search_inc
       lon_reach = lon_reach + diag_vel_search_inc
 
 
       do j=1,ny_global
       do i=1,nx_global
 
         if (   (lat_loc(nloc) - lat_reach <= ULAT_G(i,j) )   .and.  &
                (lat_loc(nloc) + lat_reach >= ULAT_G(i,j) )   .and.  &
                (lon_loc(nloc) - lon_reach <= ULON_G(i,j) )   .and.  &
                (lon_loc(nloc) + lon_reach >= ULON_G(i,j) ) ) then
 
               !-----------------------------------------------------
               ! accept point only if it is an ocean point
               !-----------------------------------------------------
 
               !---------------------------------------------------------
               !*** note that REGION_MASK_G is on the t-grid; should test
               !***  on the u-grid
               !---------------------------------------------------------

               if (REGION_MASK_G(i,j) /= 0) then

                  !---------------------------------------------------------
                  !*** note that this is crude approximation; will be changed
                  !***  in a more permanent version at a later time
                  distance = sqrt( (lat_loc(nloc)-ULAT_G(i,j))**2 + &
                                   (lon_loc(nloc)-ULON_G(i,j))**2 )
                  !---------------------------------------------------------

                  if (distance < min_distance) then
                    j_loc(nloc) = j
                    i_loc(nloc) = i
                    min_distance = distance
                    vlat_loc(nloc) = ULAT_G(i,j)
                    vlon_loc(nloc) = ULON_G(i,j)
                  endif

               endif ! ocean point
         endif

       enddo ! i
       enddo ! j

     enddo reach_loop

   enddo ! nloc

   write(stdout,*) ' '
   write(stdout,*) '   Velocity diagnostics will be computed at the following locations:'
   do nloc = 1, num_vel_loc
   write(stdout,'(2x,a,i3,a,i3,a, 3x, a,F25.15,a,F25.15,a)' ) &
                '(',   i_loc(nloc), ',',   j_loc(nloc),')',   &
                '(',vlat_loc(nloc), ',',vlon_loc(nloc),')'
   enddo
   call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)

   endif
 
   call broadcast_array(j_loc , master_task)
   call broadcast_array(i_loc , master_task)

   do nloc = 1, num_vel_loc
     if ( j_loc(nloc) == search_error .or. i_loc(nloc) == search_error)  &
       call exit_POP (sigAbort, 'velocity diagnostics lat/lon search failed')
   enddo

   deallocate (ULAT_G, ULON_G)
   deallocate (REGION_MASK_G)
  endif !ldiag_velocity

   if (my_task == master_task) then
      write(stdout,blank_fmt)
      write(stdout,delim_fmt)
      write(stdout,blank_fmt)
      write(stdout,'(a18)') 'End Diagnostic options'
      write(stdout,blank_fmt)
      write(stdout,delim_fmt)
      call POP_IOUnitsFlush(POP_stdout);  call POP_IOUnitsFlush(stdout) ! temporary
   endif



!EOC

 end subroutine init_diagnostics

!***********************************************************************
!BOP
! !IROUTINE: diag_init_sums
! !INTERFACE:

 subroutine diag_init_sums

! !DESCRIPTION:
!  Sets diagnostic flags and initializes diagnostic sums at beginning
!  of each step.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  output formats
!
!-----------------------------------------------------------------------

   character (8), parameter :: &
      time_fmt1 = '(a15,i8)'

   character (11), parameter :: &
      time_fmt2 = '(a15,f14.6)'

   character (20), parameter :: &
      date_fmt1 = '(a15,a2,1x,a3,1x,a4)'

!-----------------------------------------------------------------------
!
!  set logical switches for diagnostic timesteps
!  do not compute on first pass of Matsuno
!
!-----------------------------------------------------------------------

   ldiag_global    = check_time_flag(diag_global_flag)
   ldiag_cfl       = check_time_flag(diag_cfl_flag)
   ldiag_transport = check_time_flag(diag_transp_flag)

   if (mix_pass == 1) then
      ldiag_global    = .false.
      ldiag_cfl       = .false.
      ldiag_transport = .false.
   endif

!-----------------------------------------------------------------------
!
!  initialize arrays for global diagnostic sums.
!
!-----------------------------------------------------------------------

   if (ldiag_global) then
      DIAG_KE_ADV_2D  = c0
      DIAG_KE_HMIX_2D = c0
      DIAG_KE_VMIX_2D = c0
      DIAG_KE_PRESS_2D = c0
      DIAG_PE_2D = c0
      DIAG_TRACER_ADV_2D = c0
      DIAG_TRACER_HDIFF_2D = c0
      DIAG_TRACER_VDIFF_2D = c0
      DIAG_TRACER_SOURCE_2D = c0
      DIAG_TRACER_SFC_FLX   = c0
   endif

!-----------------------------------------------------------------------
!
!  write some stuff to log file (stdout)
!
!-----------------------------------------------------------------------

   if ((ldiag_global .or. ldiag_cfl) .and. &
       my_task == master_task) then
      write(stdout,blank_fmt)
      write(stdout,time_fmt1) 'Step number  : ',nsteps_total
      write(stdout,date_fmt1) 'Date         : ',cday,cmonth3,cyear
      write(stdout,time_fmt1) 'Hour         : ',ihour
      write(stdout,time_fmt1) 'Minute       : ',iminute
      write(stdout,time_fmt1) 'Second       : ',isecond

      if (tsecond < seconds_in_hour) then
         write(stdout,time_fmt2) 'Time(seconds): ', tsecond
      elseif (tsecond < seconds_in_day) then
         write(stdout,time_fmt2) 'Time(hours)  : ', &
                                             tsecond/seconds_in_hour
      elseif (tsecond < 31536000._r8) then
         write(stdout,time_fmt2) 'Time(days)   : ', tday
      else
         write(stdout,time_fmt2) 'Time(years)  : ', tyear
      endif
   endif

!-----------------------------------------------------------------------
!EOC

 end subroutine diag_init_sums

!***********************************************************************
!BOP
! !IROUTINE: diag_global_preupdate
! !INTERFACE:

 subroutine diag_global_preupdate(DH,DHU)

! !DESCRIPTION:
!  Computes diagnostic sums for those diagnostics computed from
!  current time prognostic variables (before updating to new time)
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
      intent(in) :: &
      DH, DHU            ! change in surface height (-Fw) at T,U pts

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      i, k,              & ! horiz, vertical level index
      n ,                & ! tracer index
      iblock ,           & ! block index
      dcount ,           & ! diag counter
      ib,ie,jb,je

   real (r8), dimension(max_blocks_clinic,6*nt+9) :: &
      local_sums           ! array for holding block sums
                           ! of each diagnostic

   real (r8) ::            &
      sum_tmp              ! temp for local sum

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1,              &! local work space
      TFACT,              &! factor for normalizing sums
      UFACT                ! factor for normalizing sums

   type (block) ::         &
      this_block           ! block information for current block

!-----------------------------------------------------------------------
!
!  compute these diagnostics only if it is time
!
!-----------------------------------------------------------------------

   if (ldiag_global) then

!-----------------------------------------------------------------------
!
!     compute local sums of each diagnostic
!     do all block sums first for better OpenMP performance
!
!-----------------------------------------------------------------------

      local_sums = c0

      !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je,i, &
      !$OMP                     k,n,dcount,WORK1,UFACT,TFACT)

      do iblock = 1,nblocks_clinic

         this_block = get_block(blocks_clinic(iblock),iblock)  

         ib = this_block%ib
         ie = this_block%ie
         jb = this_block%jb
         je = this_block%je

         TFACT = TAREA(:,:,iblock)*RCALCT(:,:,iblock)
         UFACT = UAREA(:,:,iblock)*RCALCU(:,:,iblock)
         if (ltripole_grid .and. this_block%jblock == nblocks_y) then
            do i=1,nx_block
               if (this_block%i_glob(i) > nx_global/2) UFACT(i,je) = c0
            end do
         endif

         dcount = 0

         dcount = dcount + 1
         WORK1 = (UVEL(:,:,1,curtime,iblock)*SMF(:,:,1,iblock)   &
                 +VVEL(:,:,1,curtime,iblock)*SMF(:,:,2,iblock))*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))
 
         dcount = dcount + 1
         WORK1 = DIAG_KE_ADV_2D(:,:,iblock)*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = DIAG_KE_HMIX_2D(:,:,iblock)*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = DIAG_KE_VMIX_2D(:,:,iblock)*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = DIAG_KE_PRESS_2D(:,:,iblock)*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = DIAG_PE_2D(:,:,iblock)*TFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = -(UBTROP(:,:,curtime,iblock)*  &
                   GRADPX(:,:,curtime,iblock) + &
                   VBTROP(:,:,curtime,iblock)*  &
                   GRADPY(:,:,curtime,iblock))* &
                   HU(:,:,iblock)*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = p5*(UVEL(:,:,1,curtime,iblock)**2 + &
                     VVEL(:,:,1,curtime,iblock)**2)* &
                 DHU(:,:,iblock)*UFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         dcount = dcount + 1
         WORK1 = PSURF(:,:,curtime,iblock)*DH(:,:,iblock)*TFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         do n=1,nt
            dcount = dcount + 1
            WORK1 = DIAG_TRACER_ADV_2D(:,:,n,iblock)*TFACT
            local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

            dcount = dcount + 1
            WORK1 = DIAG_TRACER_HDIFF_2D(:,:,n,iblock)*TFACT
            local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

            dcount = dcount + 1
            WORK1 = DIAG_TRACER_VDIFF_2D(:,:,n,iblock)*TFACT
            local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

            dcount = dcount + 1
            WORK1 = DIAG_TRACER_SOURCE_2D(:,:,n,iblock)*TFACT
            local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

            dcount = dcount + 1
            do k=1,km
               if (sfc_layer_type == sfc_layer_varthick .and. &
                   k == 1) then
                  WORK1 = merge(TAREA(:,:,iblock)*(                    &
                            (dz(1) + PSURF(:,:,newtime,iblock)/grav)*  &
                                TRACER(:,:,k,n,newtime,iblock) -       &
                            (dz(1) + PSURF(:,:,mixtime,iblock)/grav)*  &
                                TRACER(:,:,k,n,oldtime,iblock)),       &
                                c0, k <= KMT(:,:,iblock))
               else
                  if (partial_bottom_cells) then
                     WORK1 = merge(DZT(:,:,k,iblock)*TAREA(:,:,iblock)*&
                                   (TRACER(:,:,k,n,newtime,iblock) -   &
                                    TRACER(:,:,k,n,oldtime,iblock)),   &
                                   c0, k <= KMT(:,:,iblock))
                  else
                     WORK1 = merge(dz(k)*TAREA(:,:,iblock)*            &
                                   (TRACER(:,:,k,n,newtime,iblock) -   &
                                    TRACER(:,:,k,n,oldtime,iblock)),   &
                                   c0, k <= KMT(:,:,iblock))
                  endif
               endif

               local_sums(iblock,dcount  ) = &
               local_sums(iblock,dcount  ) + &
                                      sum(WORK1(ib:ie,jb:je))/c2dtt(k)
               local_sums(iblock,dcount+1) = &
               local_sums(iblock,dcount+1) + &
                                 sum(abs(WORK1(ib:ie,jb:je)))/c2dtt(k)
            end do  !k loop
            dcount = dcount + 1
         end do !tracer loop
      end do ! block loop

      !$OMP END PARALLEL DO 

!-----------------------------------------------------------------------
!
!     finish up scalar diagnostics
!
!-----------------------------------------------------------------------

      dcount = 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ws         = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke_adv     = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke_hmix    = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke_vmix    = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke_press   = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_pe         = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_press_free = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke_free    = global_sum(sum_tmp,distrb_clinic)/volume_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke_psfc    = global_sum(sum_tmp,distrb_clinic)/volume_t
      do n=1,nt
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        diag_tracer_adv(n)   = &
                         global_sum(sum_tmp,distrb_clinic)/volume_t
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        diag_tracer_hdiff(n) = &
                         global_sum(sum_tmp,distrb_clinic)/volume_t
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        diag_tracer_vdiff(n) = & 
                         global_sum(sum_tmp,distrb_clinic)/volume_t
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        diag_tracer_source(n) = & 
                         global_sum(sum_tmp,distrb_clinic)/volume_t
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        dtracer_avg(n) = global_sum(sum_tmp,distrb_clinic)/volume_t
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        dtracer_abs(n) = global_sum(sum_tmp,distrb_clinic)/volume_t
      end do

!-----------------------------------------------------------------------

   endif ! ldiag_global

!-----------------------------------------------------------------------
        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(HMXL(:,:,iblock), tavg_HMXL,   iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(HMXL(:,:,iblock), tavg_HMXL_2, iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(HMXL(:,:,iblock), tavg_XMXL,   iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(HMXL(:,:,iblock), tavg_XMXL_2, iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(HMXL(:,:,iblock), tavg_TMXL, iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(KPP_HBLT(:,:,iblock), tavg_HBLT, iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(KPP_HBLT(:,:,iblock), tavg_XBLT, iblock, 1)
        end do
        !$OMP END PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(iblock)
        do iblock=1,nblocks_clinic
          call accumulate_tavg_field(KPP_HBLT(:,:,iblock), tavg_TBLT, iblock, 1)
        end do
        !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!
!     accumulate movie diagnostics if requested
!
!-----------------------------------------------------------------------

      if (movie_requested(movie_HMXL) ) then
        !$OMP PARALLEL DO
        do iblock=1,nblocks_clinic
          call update_movie_field(HMXL(:,:,iblock), movie_HMXL, iblock, 1)
        end do
        !$OMP END PARALLEL DO
      endif

      if (movie_requested(movie_HBLT) ) then
        !$OMP PARALLEL DO
        do iblock=1,nblocks_clinic
          call update_movie_field(KPP_HBLT(:,:,iblock), movie_HBLT, iblock, 1)
        end do
        !$OMP END PARALLEL DO
      endif

!EOC

 end subroutine diag_global_preupdate

!***********************************************************************
!BOP
! !IROUTINE: diag_global_afterupdate
! !INTERFACE:

 subroutine diag_global_afterupdate

! !DESCRIPTION:
!  Computes diagnostic sums for those diagnostics computed from
!  updated prognostic variables (after updating to new time).
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (POP_i4) :: &
      errorCode         ! returned error code

   integer (int_kind) :: &
      i, k,              & ! horiz, vertical level index
      n ,                & ! tracer index
      iblock ,           & ! block index
      dcount ,           & ! diag counter
      ib,ie,jb,je

   real (r8), dimension(max_blocks_clinic,2*nt+2) :: &
      local_sums           ! array for holding block sums
                           ! of each diagnostic

   real (r8), dimension(max_blocks_clinic,km,nt) :: &
      local_sums_k         ! array for holding block sums
                           ! of each diagnostic at each level

   real (r8) ::            &
      sum_tmp              ! temp for local sum

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1,              &! local work space
      TFACT,              &! factor for normalizing sums
      UFACT                ! factor for normalizing sums

   type (block) ::         &
      this_block           ! block information for current block

!-----------------------------------------------------------------------
!
!  compute these diagnostics only if it is time
!
!-----------------------------------------------------------------------

   if (ldiag_global) then

!-----------------------------------------------------------------------
!
!     compute local block sums of each diagnostic
!
!-----------------------------------------------------------------------

      local_sums = c0
      if (diag_all_levels) local_sums_k = c0

      !$OMP PARALLEL DO &
      !$OMP PRIVATE(iblock,this_block,ib,ie,jb,je,i,k,n,dcount, &
      !$OMP         WORK1,TFACT,UFACT)

      do iblock = 1,nblocks_clinic

         this_block = get_block(blocks_clinic(iblock),iblock)  

         ib = this_block%ib
         ie = this_block%ie
         jb = this_block%jb
         je = this_block%je

         TFACT = TAREA(:,:,iblock)*RCALCT(:,:,iblock)
         UFACT = UAREA(:,:,iblock)*RCALCU(:,:,iblock)
         if (ltripole_grid .and. this_block%jblock == nblocks_y) then
            do i=1,nx_block
               if (this_block%i_glob(i) > nx_global/2) UFACT(i,je) = c0
            end do
         endif

         dcount = 0

         !*** global mean sea level

         dcount = dcount + 1
         WORK1 = PSURF(:,:,curtime,iblock)/grav*TFACT
         local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         !*** avg horiz KE at new time

         dcount = dcount + 1
         do k = 1,km
            if (partial_bottom_cells) then
               WORK1 = merge(p5*(UVEL(:,:,k,curtime,iblock)**2 +     &
                                 VVEL(:,:,k,curtime,iblock)**2)*     &
                             UFACT*DZU(:,:,k,iblock),                &
                             c0, KMU(:,:,iblock) >= k)
            else
               WORK1 = merge(p5*(UVEL(:,:,k,curtime,iblock)**2 +     &
                                 VVEL(:,:,k,curtime,iblock)**2)*     &
                             UFACT*dz(k),                            &
                             c0, KMU(:,:,iblock) >= k)
            endif

            local_sums(iblock,dcount) = local_sums(iblock,dcount) + & 
                                        sum(WORK1(ib:ie,jb:je))

         enddo ! k loop

         !*** average tracers at new time

         do n=1,nt
            dcount = dcount + 1
            do k = 1,km
               if (sfc_layer_type == sfc_layer_varthick .and. &
                   k == 1) then
                  WORK1 = merge(TRACER(:,:,k,n,curtime,iblock)*        &
                                TAREA(:,:,iblock)*                     &
                                (c1 + PSURF(:,:,curtime,iblock)/       &
                                (dz(1)*grav))*dz(1)                    &
                               ,c0,KMT(:,:,iblock) >= k)
               else
                  if (partial_bottom_cells) then
                     WORK1 = merge(TRACER(:,:,k,n,curtime,iblock)*     &
                                   TAREA(:,:,iblock)*DZT(:,:,k,iblock),&
                                   c0,KMT(:,:,iblock) >= k)
                  else
                     WORK1 = merge(TRACER(:,:,k,n,curtime,iblock)*  &
                                   TAREA(:,:,iblock)*dz(k),         &
                                   c0,KMT(:,:,iblock) >= k)
                  endif
               endif

               local_sums(iblock,dcount) = local_sums(iblock,dcount) + & 
                                           sum(WORK1(ib:ie,jb:je))

               if (diag_all_levels) then
                  local_sums_k(iblock,k,n) = sum(WORK1(ib:ie,jb:je))
               endif
            end do ! k loop
         end do ! tracer loop

         !*** surface tracer fluxes

         do n=1,nt
            dcount = dcount + 1

            WORK1 = DIAG_TRACER_SFC_FLX(:,:,n,iblock)*TFACT

            local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je))

         end do ! tracer loop
      end do ! block loop

      !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!
!     finish up by computing global sums of local sums
!
!-----------------------------------------------------------------------

      dcount= 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_sealevel = global_sum(sum_tmp,distrb_clinic)/area_t
      dcount = dcount + 1
      sum_tmp = sum(local_sums(:,dcount))
      diag_ke       = global_sum(sum_tmp,distrb_clinic)/volume_u
      do n=1,nt
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        avg_tracer(n) = global_sum(sum_tmp,distrb_clinic)/volume_t
        if (diag_all_levels) then
           do k=1,km
              sum_tmp = sum(local_sums_k(:,k,n))
              avg_tracer_k(k,n) = global_sum(sum_tmp,distrb_clinic)/ &
                                  area_t_k(k)
           end do
        endif
      end do
      do n=1,nt
        dcount = dcount + 1
        sum_tmp = sum(local_sums(:,dcount))
        sfc_tracer_flux(n) = global_sum(sum_tmp,distrb_clinic)/area_t
      end do

!-----------------------------------------------------------------------
!
!     convert some tracer diagnostics to different units
!
!-----------------------------------------------------------------------

      avg_tracer(2) = avg_tracer(2)*salt_to_ppt ! salinity to ppt
      if (diag_all_levels) then
        do k=1,km
          avg_tracer_k(k,2) = avg_tracer_k(k,2)*salt_to_ppt
        end do
      endif

      sfc_tracer_flux(1) = sfc_tracer_flux(1)/hflux_factor ! flux in W/m2
      sfc_tracer_flux(2) = sfc_tracer_flux(2)/salinity_factor ! FW flux

!-----------------------------------------------------------------------
!
!     retrieve solver diagnostics
!
!-----------------------------------------------------------------------

      call POP_SolversGetDiagnostics(iterationCount, rmsResidual, &
                                     errorCode)

!-----------------------------------------------------------------------

   endif ! ldiag_global

!-----------------------------------------------------------------------
!EOC

 end subroutine diag_global_afterupdate

!***********************************************************************
!BOP
! !IROUTINE: diag_print
! !INTERFACE:

 subroutine diag_print

! !DESCRIPTION:
!  Writes scalar diagnostic info to both stdout and diagnostic output
!  file.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variable
!
!-----------------------------------------------------------------------

   integer (int_kind) ::  &
      k,n,                &! dummy indices
      ier

   character (25) :: diag_name  ! name of output diagnostic

   character (29), parameter :: &
      diag_fmt = '(1pe22.15,1x,1pe22.15,2x,a25)'

   character (14), parameter :: &
      stdo_fmt = '(a32,1pe22.15)'

   character (21), parameter :: &
      trcr_fmt = '(a24,i3,a19,1pe22.15)'

   character (char_len_long) ::  &
      string

!-----------------------------------------------------------------------
!
!  print diagnostics if it is time
!
!-----------------------------------------------------------------------

   if (ldiag_global) then
      if (my_task == master_task) then

         !*** append diagnostic output to end of diag output file

         open(diag_unit, file=diag_outfile, status='old', &
                         position='append')

!-----------------------------------------------------------------------
!
!        print global diagnostics to stdout
!
!-----------------------------------------------------------------------

         write(stdout,blank_fmt)
         write(stdout,'(a17)') 'K.E. diagnostics:'
         write(stdout,stdo_fmt) &
         ' mean K.E.                     :', diag_ke
         write(stdout,stdo_fmt) &
         ' K.E. change from horiz. visc. :', diag_ke_hmix
         write(stdout,stdo_fmt) &
         ' K.E. change from vert.  visc. :', diag_ke_vmix
         write(stdout,stdo_fmt) &
         ' K.E. change from wind stress  :', diag_ws

         write(stdout,blank_fmt)
         write(stdout,stdo_fmt) &
         ' advective K.E. balance        :', diag_ke_adv
         write(stdout,stdo_fmt) &
         '                               :', diag_ke_free
         write(stdout,stdo_fmt) &
         ' P.E. verses pressure work     :', diag_pe
         write(stdout,stdo_fmt) &
         '                               :', diag_ke_press
         write(stdout,stdo_fmt) &
         ' surface-pressure contribution :', diag_ke_psfc
         write(stdout,stdo_fmt) &
         '                               :', diag_press_free

         write(stdout,blank_fmt)
         write(stdout,*) 'Tracer diagnostics:'
         do n=1,nt
            write(stdout,blank_fmt)
            write(stdout,trcr_fmt) &
            ' mean value  of tracer  ',n, &
            '                  : ',    avg_tracer(n)
            write(stdout,trcr_fmt) &
            ' abs  change in tracer  ',n, &
            '                  : ',    dtracer_abs(n)
            write(stdout,trcr_fmt) &
            ' mean change in tracer  ',n, &
            '                  : ',    dtracer_avg(n)
            write(stdout,trcr_fmt) &
            ' mean change in tracer  ',n, &
            ' from advection   : ',    diag_tracer_adv(n)
            write(stdout,trcr_fmt) &
            ' mean change in tracer  ',n, &
            ' from horiz. diff.: ',    diag_tracer_hdiff(n)
            write(stdout,trcr_fmt) &
            ' mean change in tracer  ',n, &
            ' from vert.  diff.: ',    diag_tracer_vdiff(n)
            write(stdout,trcr_fmt) &
            ' mean change in tracer  ',n, &
            ' from sources     : ',    diag_tracer_source(n)
            write(stdout,trcr_fmt) &
            ' surface flux of tracer ',n, &
            '                  : ',    sfc_tracer_flux(n)

            if (diag_all_levels) then
               do k=1,km
                   write(stdout,'(a22,i3,a10,i3,a3,1pe22.15)') &
                         ' mean value of tracer ', n,          &
                         ' at level ',             k,          &
                         ' : ',                    avg_tracer_k(k,n)
               end do
            endif
         end do

         write(stdout,blank_fmt)
         write(stdout,'(a18)') 'Other diagnostics:'
         write(stdout,'(a24,1pe22.15)') ' global mean sea level: ', &
                                          diag_sealevel
         write(stdout,'(a19,1pe22.15,2x,i4)') ' residual, scans : ', &
                                          rmsResidual, iterationCount

!-----------------------------------------------------------------------
!
!        print global diagnostics to output file
!
!-----------------------------------------------------------------------

         diag_name = 'mean KE'
         write(diag_unit,diag_fmt) tday, diag_ke, diag_name

         do n=1,nt
            write(diag_name,'(a12,i2)') 'mean tracer ',n
            write(diag_unit,diag_fmt) tday, avg_tracer(n), diag_name
            if (diag_all_levels) then
               do k=1,km
                  write(diag_name,'(a12,i2,a8,i2)') &
                        'mean tracer ',n,' at lvl ',k
                  write(diag_unit,diag_fmt) tday, avg_tracer_k(k,n), &
                                            diag_name
               end do
            endif
         end do

         diag_name = 'mean sea level'
         write(diag_unit,diag_fmt) tday, diag_sealevel, diag_name
         diag_name = 'solver residual'
         write(diag_unit,diag_fmt) tday, rmsResidual, diag_name
         diag_name = 'solver iterations'
         write(diag_unit,diag_fmt) tday, real(iterationCount), diag_name

         do n=1,nt
            write(diag_name,'(a7,i3,a11)') 'Tracer ',n,' change avg'
            write(diag_unit,diag_fmt) tday, dtracer_avg(n), diag_name
            write(diag_name,'(a7,i3,a11)') 'Tracer ',n,' change abs'
            write(diag_unit,diag_fmt) tday, dtracer_abs(n), diag_name

            write(diag_name,'(a7,i3,a14)') 'Tracer ',n,' change advect'
            write(diag_unit,diag_fmt) tday, diag_tracer_adv(n), diag_name
            write(diag_name,'(a7,i3,a13)') 'Tracer ',n,' change hdiff'
            write(diag_unit,diag_fmt) tday, diag_tracer_hdiff(n), diag_name
            write(diag_name,'(a7,i3,a13)') 'Tracer ',n,' change vdiff'
            write(diag_unit,diag_fmt) tday, diag_tracer_vdiff(n), diag_name
            write(diag_name,'(a7,i3,a14)') 'Tracer ',n,' change source'
            write(diag_unit,diag_fmt) tday, diag_tracer_source(n), diag_name

            write(diag_name,'(a7,i3,a14)') 'Tracer ',n,' surface flux'
            write(diag_unit,diag_fmt) tday, sfc_tracer_flux(n), diag_name
         end do

         diag_name = 'KE change horiz visc'
         write(diag_unit,diag_fmt) tday, diag_ke_hmix, diag_name
         diag_name = 'KE change vert  visc'
         write(diag_unit,diag_fmt) tday, diag_ke_vmix, diag_name
         diag_name = 'KE change wind stress'
         write(diag_unit,diag_fmt) tday, diag_ws, diag_name

         diag_name = 'advect KE balance'
         write(diag_unit,diag_fmt) tday, diag_ke_adv, diag_name
         diag_name = 'advect KE bal free sfc'
         write(diag_unit,diag_fmt) tday, diag_ke_free, diag_name
         diag_name = 'PE change'
         write(diag_unit,diag_fmt) tday, diag_pe, diag_name
         diag_name = 'pressure work'
         write(diag_unit,diag_fmt) tday, diag_ke_press, diag_name
         diag_name = 'KE change sfc-pressure'
         write(diag_unit,diag_fmt) tday, diag_ke_psfc, diag_name
         diag_name = 'press work free sfc'
         write(diag_unit,diag_fmt) tday, diag_press_free, diag_name

!-----------------------------------------------------------------------
!
!        close file and print timer info
!
!-----------------------------------------------------------------------

         close(diag_unit)

      endif ! master_task

      call timer_print_all()
      call document ('diag_print', 'file written: '// trim(diag_outfile))


   endif ! ldiag_global


!-----------------------------------------------------------------------
!EOC

 end subroutine diag_print

!***********************************************************************
!BOP
! !IROUTINE: diag_transport
! !INTERFACE:

 subroutine diag_transport

! !DESCRIPTION:
!  Calculates and prints various user-defined transports.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      i,j,k,n,iblock,    &! dummy loop indices
      ib,ie,jb,je,       &! beg,end indices for block physical domain
      ier

   real (r8), dimension(nx_block,ny_block) :: &
      MASS_Z,              &! zonal transport of mass
      HEAT_Z,              &! zonal transport of heat
      SALT_Z,              &! zonal transport of salt
      MASS_M,              &! merid transport of mass
      HEAT_M,              &! merid transport of heat
      SALT_M,              &! merid transport of salt
      WORK1, WORK2          ! local temp space

   real (r8) ::            &
      sum_tmp               ! temp for computing block sums

   real (r8), dimension(:,:), allocatable ::            &
      mass_tran,           &! local sum of mass transport
      heat_tran,           &! local sum of heat transport
      salt_tran             ! local sum of salt transport

   character (char_len_long) ::  &
      string

   type (block) ::         &
      this_block           ! block information for current block

!-----------------------------------------------------------------------
!
!  only compute if it is time
!  initialize transport sums
!
!-----------------------------------------------------------------------

   if (ldiag_transport) then

      if (my_task == master_task) then
         !*** re-open file and append current values to previous
         !*** outputs from this run

         open(trans_unit, file=diag_transport_outfile, &
              status='old', position='append')
      endif

!-----------------------------------------------------------------------
!
!     calculate transports for each block
!
!-----------------------------------------------------------------------

      allocate(mass_tran(nblocks_clinic,num_transports), &
               heat_tran(nblocks_clinic,num_transports), &
               salt_tran(nblocks_clinic,num_transports))

      !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je,k,j,i,n, &
      !$OMP                     WORK1,WORK2,MASS_M,SALT_M,HEAT_M, &
      !$OMP                     MASS_Z,SALT_Z,HEAT_Z)

      do iblock = 1,nblocks_clinic

         this_block = get_block(blocks_clinic(iblock),iblock)  

         ib = this_block%ib
         ie = this_block%ie
         jb = this_block%jb
         je = this_block%je

         mass_tran(iblock,:) = c0
         heat_tran(iblock,:) = c0
         salt_tran(iblock,:) = c0

         do k = 1,km

            if (partial_bottom_cells) then
               do j=jb,je
               do i=ib,ie
                  MASS_M(i,j) = p5*(                          &
                                UVEL(i,j  ,k,curtime,iblock)* &
                                DYU(i,j  ,iblock)*            &
                                DZU(i,j  ,k,iblock) +         &
                                UVEL(i,j-1,k,curtime,iblock)* &
                                DYU(i,j-1,iblock)*            &
                                DZU(i,j-1,k,iblock) )

                  MASS_Z(i,j) = p5*(                             &
                                VVEL(i  ,j,k,curtime,iblock)*    &
                                DXU(i  ,j,iblock)*               &
                                DZU(i  ,j,k,iblock) +            &
                                VVEL(i-1,j,k,curtime,iblock)*    &
                                DXU(i-1,j,iblock)*               &
                                DZU(i-1,j,k,iblock) )

               end do
               end do
            else
               do j=jb,je
               do i=ib,ie
                  MASS_M(i,j) = p5*(                                   &
                                UVEL(i,j  ,k,curtime,iblock)*          &
                                DYU(i,j  ,iblock)*dz(k)                &
                              + UVEL(i,j-1,k,curtime,iblock)*          &
                                DYU(i,j-1,iblock)*dz(k) )

                  MASS_Z(i,j) = p5*(                                   &
                                VVEL(i  ,j,k,curtime,iblock)*          &
                                DXU(i  ,j,iblock)*dz(k)                &
                              + VVEL(i-1,j,k,curtime,iblock)*          &
                                DXU(i-1,j,iblock)*dz(k) )
               end do
               end do
            endif

            do j=jb,je
            do i=ib,ie
               HEAT_M(i,j) = p5*MASS_M(i,j)*                     &
                             (TRACER(i+1,j,k,1,curtime,iblock) + &
                              TRACER(i  ,j,k,1,curtime,iblock))

               SALT_M(i,j) = p5*MASS_M(i,j)*                     &
                             (TRACER(i+1,j,k,2,curtime,iblock) + &
                              TRACER(i  ,j,k,2,curtime,iblock))

               HEAT_Z(i,j) = p5*MASS_Z(i,j)*                     &
                             (TRACER(i,j+1,k,1,curtime,iblock) + &
                              TRACER(i,j  ,k,1,curtime,iblock))

               SALT_Z(i,j) = p5*MASS_Z(i,j)*                     &
                             (TRACER(i,j+1,k,2,curtime,iblock) + &
                              TRACER(i,j  ,k,2,curtime,iblock))
            end do
            end do

!-----------------------------------------------------------------------
!
!           compute transport sums for each defined transport
!
!-----------------------------------------------------------------------

            do n=1,num_transports
               if (k >= transports(n)%add(5,iblock) .and.  &
                   k <= transports(n)%add(6,iblock) ) then

                  !*** local sum

                  select case (transports(n)%type)

                  case (transport_type_zonal)

                     do j=transports(n)%add(3,iblock), &
                          transports(n)%add(4,iblock)
                     do i=transports(n)%add(1,iblock), &
                          transports(n)%add(2,iblock)
                        mass_tran(iblock,n) = mass_tran(iblock,n) + &
                                              MASS_Z(i,j)
                        heat_tran(iblock,n) = heat_tran(iblock,n) + &
                                              HEAT_Z(i,j)
                        salt_tran(iblock,n) = salt_tran(iblock,n) + &
                                              SALT_Z(i,j)
                     end do
                     end do

                  case (transport_type_merid)

                     do j=transports(n)%add(3,iblock), &
                          transports(n)%add(4,iblock)
                     do i=transports(n)%add(1,iblock), &
                          transports(n)%add(2,iblock)
                        mass_tran(iblock,n) = mass_tran(iblock,n) + &
                                              MASS_M(i,j)
                        heat_tran(iblock,n) = heat_tran(iblock,n) + &
                                              HEAT_M(i,j)
                        salt_tran(iblock,n) = salt_tran(iblock,n) + &
                                              SALT_M(i,j)
                     end do
                     end do

                  end select

               endif
            end do  ! transport loop
         end do ! k loop
      end do ! block loop

      !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!
!     finish global sums for each transport, convert to useful units
!     and write to output file
!
!-----------------------------------------------------------------------

      do n=1,num_transports

         sum_tmp = sum(mass_tran(:,n))
         transports(n)%mass = global_sum(sum_tmp,distrb_clinic)* &
                              mass_to_Sv

         sum_tmp = sum(heat_tran(:,n))
         transports(n)%heat = global_sum(sum_tmp,distrb_clinic)* &
                              heat_to_PW

         sum_tmp = sum(salt_tran(:,n))
         transports(n)%salt = global_sum(sum_tmp,distrb_clinic)* &
                              salt_to_Svppt

         if (my_task == master_task) then
            write(trans_unit,'(4(1pe13.6,2x), a)')           &
                                   tday, transports(n)%mass, &
                                         transports(n)%heat, &
                                         transports(n)%salt, &
                                         trim(transports(n)%name)
         endif
      end do ! transport loop

      if (my_task == master_task) then
        close(trans_unit)
      endif ! master_task

      call document ('diag_transport', 'file written: '//trim(diag_transport_outfile))

      deallocate(mass_tran, heat_tran, salt_tran)

   endif ! ldiag_transport

!-----------------------------------------------------------------------
!EOC

 end subroutine diag_transport

!***********************************************************************
!BOP
! !IROUTINE: cfl_advect
! !INTERFACE:

 subroutine cfl_advect(k,iblock,UTE,UTW,VTN,VTS,WTK,WTKB,this_block)

! !DESCRIPTION:
!  Calculates maximum advective cfl number for a particular block.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) :: &
      k,          &! vertical level index
      iblock       ! local block index in baroclinic distribution

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      UTE, UTW,          &! advective fluxes through east, west  faces
      VTN, VTS,          &! advective fluxes through north,south faces
      WTK, WTKB           ! advective fluxes through top, bottom faces

   type (block), intent(in) ::     &
      this_block           ! block information for current block

!EOP
!BOC
!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      i,j                ! dummy loop indices

   real (r8), dimension(nx_block,ny_block) :: &
      CFL_E, CFL_W, &! cfl numbers due to advective fluxes through
      CFL_N, CFL_S, &! each face
      CFL_T, CFL_B, &! 
      CFL_TOT        ! sum for total cfl number

!-----------------------------------------------------------------------
!
!  only compute cfl numbers if it is time
!
!-----------------------------------------------------------------------

   if (ldiag_cfl) then

!-----------------------------------------------------------------------
!
!     advective fluxes through east,west,north,south,top,bottom faces
!
!-----------------------------------------------------------------------

      if (partial_bottom_cells) then

         CFL_E = abs(UTE*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k)
         CFL_W = abs(UTW*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k)
         CFL_N = abs(VTN*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k)
         CFL_S = abs(VTS*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k)
         CFL_T = abs(WTK /DZT(:,:,k,iblock))*dt(k)
         CFL_B = abs(WTKB/DZT(:,:,k,iblock))*dt(k)

      else

         CFL_E = abs(UTE*TAREA_R(:,:,iblock))*dt(k)
         CFL_W = abs(UTW*TAREA_R(:,:,iblock))*dt(k)
         CFL_N = abs(VTN*TAREA_R(:,:,iblock))*dt(k)
         CFL_S = abs(VTS*TAREA_R(:,:,iblock))*dt(k)
         CFL_T = abs(WTK /dz(k))*dt(k)
         CFL_B = abs(WTKB/dz(k))*dt(k)

      endif

      CFL_TOT = p5*(CFL_E + CFL_W + CFL_N + CFL_S + CFL_T + CFL_B)

!-----------------------------------------------------------------------
!
!     find max cfl numbers and locations for this block
!
!-----------------------------------------------------------------------

      cfl_advuk_block(iblock,k) = -c1
      cfl_advvk_block(iblock,k) = -c1
      cfl_advwk_block(iblock,k) = -c1
      cfl_advtk_block(iblock,k) = -c1
      cfladd_advuk_block(:,iblock,k) = 0
      cfladd_advvk_block(:,iblock,k) = 0
      cfladd_advwk_block(:,iblock,k) = 0
      cfladd_advtk_block(:,iblock,k) = 0

      do j=this_block%jb,this_block%je
      do i=this_block%ib,this_block%ie

         !*** zonal cfl number east flux

         if (CFL_E(i,j) > cfl_advuk_block(iblock,k)) then
            cfl_advuk_block(iblock,k) = CFL_E(i,j)
            cfladd_advuk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advuk_block(2,iblock,k) = this_block%j_glob(j)
         endif

         !*** zonal cfl number west flux

         if (CFL_W(i,j) > cfl_advuk_block(iblock,k)) then
            cfl_advuk_block(iblock,k) = CFL_W(i,j)
            cfladd_advuk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advuk_block(2,iblock,k) = this_block%j_glob(j)
         endif

         !*** meridional cfl number north flux

         if (CFL_N(i,j) > cfl_advvk_block(iblock,k)) then
            cfl_advvk_block(iblock,k) = CFL_N(i,j)
            cfladd_advvk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advvk_block(2,iblock,k) = this_block%j_glob(j)
         endif

         !*** meridional cfl number south flux

         if (CFL_S(i,j) > cfl_advvk_block(iblock,k)) then
            cfl_advvk_block(iblock,k) = CFL_S(i,j)
            cfladd_advvk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advvk_block(2,iblock,k) = this_block%j_glob(j)
         endif

         !*** vertical cfl number top flux

         if (CFL_T(i,j) > cfl_advwk_block(iblock,k)) then
            cfl_advwk_block(iblock,k) = CFL_T(i,j)
            cfladd_advwk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advwk_block(2,iblock,k) = this_block%j_glob(j)
         endif

         !*** vertical cfl number bottom flux

         if (CFL_B(i,j) > cfl_advwk_block(iblock,k)) then
            cfl_advwk_block(iblock,k) = CFL_B(i,j)
            cfladd_advwk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advwk_block(2,iblock,k) = this_block%j_glob(j)
         endif

         !*** total advective CFL number

         if (CFL_TOT(i,j) > cfl_advtk_block(iblock,k)) then
            cfl_advtk_block(iblock,k) = CFL_TOT(i,j)
            cfladd_advtk_block(1,iblock,k) = this_block%i_glob(i)
            cfladd_advtk_block(2,iblock,k) = this_block%j_glob(j)
         endif

      end do
      end do

!-----------------------------------------------------------------------
!
!     finished computing advective cfl numbers
!
!-----------------------------------------------------------------------

   endif ! ldiag_cfl

!-----------------------------------------------------------------------
!EOC

 end subroutine cfl_advect

!***********************************************************************
!BOP
! !IROUTINE: cfl_vdiff
! !INTERFACE:

 subroutine cfl_vdiff(k,VDC,VVC,this_block)

! !DESCRIPTION:
!  Calculates maximum vertical diffusion cfl numbers and their
!  locations in i,j,k directions.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) :: &
      k              ! vertical level index

   real (r8), dimension(:,:,:,:), intent(in) :: &
      VDC                 ! tracer diffusivity

   real (r8), dimension(:,:,:), intent(in) :: &
      VVC                 ! viscosity

   type (block), intent(in) ::     &
      this_block           ! block information for current block

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) ::   &
      kmax, klvl, km1,     &! vertical level indices
      i,j,n,               &! dummy index
      bid                   ! local block id

   real (r8) :: &
      local_max             ! temp for holding local block max

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), save :: &
      VDCM1, VVCM1        ! save k-1 values for next k level

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1,             &! local work space
      DZ1, DZ2, DZ3       ! grid arrays

!-----------------------------------------------------------------------
!
!  only compute cfl numbers if it is time
!
!-----------------------------------------------------------------------

   if (ldiag_cfl) then

!-----------------------------------------------------------------------
!
!     determine size of VDC arrays and set up some grid arrays
!
!-----------------------------------------------------------------------

      bid = this_block%local_id
      kmax = size(VDC, DIM=3)

      if (kmax > km) then  ! KPP uses VDC(0:km+1)
         klvl = k+1
         km1  = k
      else
         klvl = k
         km1  = k-1
      endif

      if (partial_bottom_cells) then
         if (k < km) then
            DZ1 = c2/(DZT(:,:,k,bid) + DZT(:,:,k+1,bid))
         else
            DZ1 = c2/DZT(:,:,k,bid)
         endif
         if (k == 1) then
            DZ2 = c2/DZT(:,:,k,bid)
         else
            DZ2 = c2/(DZT(:,:,k,bid) + DZT(:,:,k-1,bid))
         endif
         DZ3 = c1/DZT(:,:,k,bid)
      else
         DZ1 = dzwr(k  )
         DZ2 = dzwr(k-1)
         DZ3 = dzr (k  )
      endif

!-----------------------------------------------------------------------
!
!     if VDC a 2d array, compute cfl number for this level and save
!     for next level
!
!-----------------------------------------------------------------------

      if (kmax == 1) then

         if (k == 1) VDCM1(:,:,bid) = c0

         WORK1 = abs(c2*(VDC(:,:,1,1)*DZ1 +  &
                         VDCM1(:,:,bid)*DZ2)*DZ3)

         !*** swap to prepare for next k level
         VDCM1(:,:,bid) = VDC(:,:,1,1)

!-----------------------------------------------------------------------
!
!     if VDC a 3d array, use k-1 value directly from array
!
!-----------------------------------------------------------------------

      else

         if (k == 1) then   !*** k=1 case (k-1 values = 0)
            WORK1 = c0
            do n=1,size(VDC,DIM=4)
               WORK1 = max(WORK1, &
                           abs(c2*VDC(:,:,klvl,n)*dzwr(1)*dzr(1)))
            end do

         else   !*** all other k levels

            WORK1 = c0
            do n=1,size(VDC,DIM=4)
               WORK1 = max(WORK1,abs(c2*(             &
                           VDC(:,:,klvl,n)*DZ1 +      &
                           VDC(:,:,km1 ,n)*DZ2)*DZ3))
            end do

         endif ! k=1

      endif ! kmax = 1

!-----------------------------------------------------------------------
!
!     find max tracer cfl for this block
!
!-----------------------------------------------------------------------

      cfl_vdifftk_block(bid,k) = c0
      local_max = -c1

      do j=this_block%jb,this_block%je
      do i=this_block%ib,this_block%ie

         if (WORK1(i,j) > local_max) then
            local_max = WORK1(i,j)
            cfladd_vdifftk_block(1,bid,k) = this_block%i_glob(i)
            cfladd_vdifftk_block(2,bid,k) = this_block%j_glob(j)
         endif

      end do
      end do

      cfl_vdifftk_block(bid,k) = local_max*dt(k)

!-----------------------------------------------------------------------
!
!     repeat for momentum diffusion
!
!-----------------------------------------------------------------------

      kmax = size(VVC, DIM=3)

      if (partial_bottom_cells) then
         if (k < km) then
            DZ1 = c2/(DZU(:,:,k,bid) + DZU(:,:,k+1,bid))
         else
            DZ1 = c2/DZU(:,:,k,bid)
         endif
         if (k == 1) then
            DZ2 = c2/DZU(:,:,k,bid)
         else
            DZ2 = c2/(DZU(:,:,k,bid) + DZU(:,:,k-1,bid))
         endif
         DZ3 = c1/DZU(:,:,k,bid)
      else
         DZ1 = dzwr(k  )
         DZ2 = dzwr(k-1)
         DZ3 = dzr (k  )
      endif

!-----------------------------------------------------------------------
!
!     if VVC a 2d array, compute cfl number for this level and save
!     for next level
!
!-----------------------------------------------------------------------

      if (kmax == 1) then

         if (k == 1) then
            VVCM1(:,:,bid) = c0
         endif

         WORK1 = abs(c2*(VVC(:,:,1)  *DZ1 + &
                         VVCM1(:,:,bid)*DZ2)*DZ3)

         VVCM1(:,:,bid) = VVC(:,:,1)    ! swap to prepare for next k level

!-----------------------------------------------------------------------
!
!     if VVC a 3d array, use k-1 value directly from array
!
!-----------------------------------------------------------------------

      else

         if (k == 1) then  !*** k=1 case (k-1 values = 0)

            WORK1 = abs(c2*VVC(:,:,1)*dzwr(1)*dzr(1))

         else  !*** all other k levels

            WORK1 = abs(c2*(VVC(:,:,k  )*DZ1 + &
                            VVC(:,:,k-1)*DZ2)*DZ3)

         endif ! k=1

      endif ! kmax=1

!-----------------------------------------------------------------------
!
!     find max momentum cfl for this block
!
!-----------------------------------------------------------------------

      cfl_vdiffuk_block(bid,k) = c0
      local_max = -c1

      do j=this_block%jb,this_block%je
      do i=this_block%ib,this_block%ie

         if (WORK1(i,j) > local_max) then
            local_max = WORK1(i,j)
            cfladd_vdiffuk_block(1,bid,k) = this_block%i_glob(i)
            cfladd_vdiffuk_block(2,bid,k) = this_block%j_glob(j)
         endif

      end do
      end do

      cfl_vdiffuk_block(bid,k) = local_max*dtu

!-----------------------------------------------------------------------
!
!  finished computing vertical diffusion cfl numbers
!
!-----------------------------------------------------------------------

   endif ! ldiag_cfl

!-----------------------------------------------------------------------
!EOC

 end subroutine cfl_vdiff

!***********************************************************************
!BOP
! !IROUTINE: cfl_hdiff
! !INTERFACE:

 subroutine cfl_hdiff(k,iblock,HDIFFCFL,t_or_u,this_block)

! !DESCRIPTION:
!  Calculates maximum horizontal diffusion cfl numbers and their
!     locations in i,j,k directions
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) :: &
      k,          &! vertical level index
      iblock,     &! local block index in baroclinic distribution
      t_or_u       ! called from hdifft(1) or hdiffu(2)

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      HDIFFCFL            ! hdiff cfl number at grid points in block

   type (block), intent(in) ::     &
      this_block           ! block information for current block

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      i,j                ! dummy loop indices

   real (r8) :: cfltemp  ! temporary for holding local maximum

   integer (int_kind), dimension(2) :: &
      cfltemp_add        ! temporary for holding address of local max

!-----------------------------------------------------------------------
!
!  compute max cfl number for this block and level
!
!-----------------------------------------------------------------------

   cfltemp = c0
   cfltemp_add(:) = 0

   do j=this_block%jb,this_block%je
   do i=this_block%ib,this_block%ie
      if (HDIFFCFL(i,j) > cfltemp) then
         cfltemp = HDIFFCFL(i,j)
         cfltemp_add(1) = this_block%i_glob(i)
         cfltemp_add(2) = this_block%j_glob(j)
      endif
   end do
   end do

   select case (t_or_u)

   case(1)
      cfl_hdifftk_block(iblock,k) = cfltemp*dt(k)
      cfladd_hdifftk_block(:,iblock,k) = cfltemp_add(:) 

   case(2)
      cfl_hdiffuk_block(iblock,k) = cfltemp*dtu
      cfladd_hdiffuk_block(:,iblock,k) = cfltemp_add(:) 

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine cfl_hdiff

!***********************************************************************
!BOP
! !IROUTINE: diag_velocity
! !INTERFACE:

 subroutine diag_velocity

! !DESCRIPTION:
!  Writes velocity diagnostics info to velocity output file.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variable
!
!-----------------------------------------------------------------------

   character (15), parameter :: &
     diag_fmt = '(4(1pe13.6,2x))'

   character (char_len_long) ::  &
     string

   integer (int_kind) :: &
     ier, n


   real (r8), dimension(nx_global,ny_global) ::  &
     ARRAY_G


   if ( lccsm ) then

     call gather_global (ARRAY_G, VVEL(:,:,1,curtime,:), master_task, &
                         distrb_clinic)

     if ( my_task == master_task ) then

       !*** append velocity output to end of velocity output file

       open (velocity_unit, file=diag_velocity_outfile, status='old', &
            position='append')

       do n=1,num_vel_loc
          write (velocity_unit,diag_fmt) tday, vlon_loc(n), vlat_loc(n), &
                                         ARRAY_G(i_loc(n),j_loc(n))
       enddo

       close (velocity_unit)

     endif 

   endif

!-----------------------------------------------------------------------
!EOC

 end subroutine diag_velocity 

!***********************************************************************
!BOP
! !IROUTINE: cfl_check
! !INTERFACE:

 subroutine cfl_check

! !DESCRIPTION:
!  Finishes the calculation of all cfl numbers by finding the
!  global maxima and locations from the local block maxima.
!  The final results are written to stdout (or log file)
!  and diagnostic output file.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      k,iblock            ! loop indices

   integer (int_kind), dimension(1) :: &
      imax                ! max block location

   real (r8) ::         &
      cfl_temp,         &! temp for holding local max
      cfl_advu,         &! max total cfl number for zonal advection
      cfl_advv,         &! max total cfl number for merid advection
      cfl_advw,         &! max total cfl number for vert  advection
      cfl_advt,         &! max total cfl number for total advection
      cfl_vdifft,       &! max cfl num for vert tracer diffusion
      cfl_vdiffu,       &! max cfl num for vert momentum diffusion
      cfl_hdifft,       &! max cfl num for horiz tracer diffusion
      cfl_hdiffu         ! max cfl num for horiz momentum diffusion

   real (r8), dimension(km) ::         & ! max for each level
      cfl_advuk,        &! max total cfl number for zonal advection
      cfl_advvk,        &! max total cfl number for merid advection
      cfl_advwk,        &! max total cfl number for vert  advection
      cfl_advtk,        &! max total cfl number for total advection
      cfl_vdifftk,      &! max cfl num for vert tracer diffusion
      cfl_vdiffuk,      &! max cfl num for vert momentum diffusion
      cfl_hdifftk,      &! max cfl num for horiz tracer diffusion
      cfl_hdiffuk        ! max cfl num for horiz momentum diffusion

   integer (int_kind), dimension(3) :: &
      cfladd_temp,      &! temp for holding local max address
      cfladd_advu,      &! address of max zonal advection cfl number
      cfladd_advv,      &! address of max merid advection cfl number
      cfladd_advw,      &! address of max vert  advection cfl number
      cfladd_advt,      &! address of max total advection cfl number
      cfladd_vdifft,    &! address of max vertical tracer diff cfl num
      cfladd_vdiffu,    &! address of max vertical moment diff cfl num
      cfladd_hdifft,    &! address of max horiz tracer diff cfl num
      cfladd_hdiffu      ! address of max horiz momentum diff cfl num

   integer (int_kind), dimension(2,km) :: & ! addresses for each level
      cfladd_advuk,     &! address of max zonal advection cfl number
      cfladd_advvk,     &! address of max merid advection cfl number
      cfladd_advwk,     &! address of max vert  advection cfl number
      cfladd_advtk,     &! address of max total advection cfl number
      cfladd_vdifftk,   &! address of max vertical tracer diff cfl num
      cfladd_vdiffuk,   &! address of max vertical moment diff cfl num
      cfladd_hdifftk,   &! address of max horiz tracer diff cfl num
      cfladd_hdiffuk     ! address of max horiz momentum diff cfl num

   character (25) :: &
      name_temp1, name_temp2  ! temporaries for writing names

   character (*), parameter :: &
      cfl_fmt1 =  "(' max. cfl number for ',a25,': ',1pe12.4,4x,3i5)", &
      cfl_fmt2 = '(1pe22.15,1x,1pe22.15,2x,a25)'

!-----------------------------------------------------------------------
!
!  compute cfl diagnostics only if it is time
!
!-----------------------------------------------------------------------

   if (ldiag_cfl) then

!-----------------------------------------------------------------------
!
!     initialize max cfl numbers
!
!-----------------------------------------------------------------------

      cfl_advu   = -c1
      cfl_advv   = -c1
      cfl_advw   = -c1
      cfl_advt   = -c1
      cfl_vdifft = -c1
      cfl_vdiffu = -c1
      cfl_hdifft = -c1
      cfl_hdiffu = -c1

      cfl_advuk   = -c1
      cfl_advvk   = -c1
      cfl_advwk   = -c1
      cfl_advtk   = -c1
      cfl_vdifftk = -c1
      cfl_vdiffuk = -c1
      cfl_hdifftk = -c1
      cfl_hdiffuk = -c1

      cfladd_advu   = 0
      cfladd_advv   = 0
      cfladd_advw   = 0
      cfladd_advt   = 0
      cfladd_vdifft = 0
      cfladd_vdiffu = 0
      cfladd_hdifft = 0
      cfladd_hdiffu = 0

      cfladd_advuk   = 0
      cfladd_advvk   = 0
      cfladd_advwk   = 0
      cfladd_advtk   = 0
      cfladd_vdifftk = 0
      cfladd_vdiffuk = 0
      cfladd_hdifftk = 0
      cfladd_hdiffuk = 0

!-----------------------------------------------------------------------
!
!     compute level maximum cfl numbers from block results
!
!-----------------------------------------------------------------------

      do k=1,km

         !*** zonal advection

         imax = maxloc(cfl_advuk_block(:,k))
         cfl_temp = cfl_advuk_block(imax(1),k)
         cfl_advuk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_advuk(k)) then
            cfladd_temp(1:2) = cfladd_advuk_block(:,imax(1),k)
         endif
         cfladd_advuk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_advuk(2,k) = global_maxval(cfladd_temp(2))

         !*** meridional advection

         imax = maxloc(cfl_advvk_block(:,k))
         cfl_temp = cfl_advvk_block(imax(1),k)
         cfl_advvk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_advvk(k)) then
            cfladd_temp(1:2) = cfladd_advvk_block(:,imax(1),k)
         endif
         cfladd_advvk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_advvk(2,k) = global_maxval(cfladd_temp(2))

         !*** vertical advection

         imax = maxloc(cfl_advwk_block(:,k))
         cfl_temp = cfl_advwk_block(imax(1),k)
         cfl_advwk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_advwk(k)) then
            cfladd_temp(1:2) = cfladd_advwk_block(:,imax(1),k)
         endif
         cfladd_advwk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_advwk(2,k) = global_maxval(cfladd_temp(2))

         !*** total advection

         imax = maxloc(cfl_advtk_block(:,k))
         cfl_temp = cfl_advtk_block(imax(1),k)
         cfl_advtk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_advtk(k)) then
            cfladd_temp(1:2) = cfladd_advtk_block(:,imax(1),k)
         endif
         cfladd_advtk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_advtk(2,k) = global_maxval(cfladd_temp(2))

         !*** vertical tracer diffusion

         imax = maxloc(cfl_vdifftk_block(:,k))
         cfl_temp = cfl_vdifftk_block(imax(1),k)
         cfl_vdifftk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_vdifftk(k)) then
            cfladd_temp(1:2) = cfladd_vdifftk_block(:,imax(1),k)
         endif
         cfladd_vdifftk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_vdifftk(2,k) = global_maxval(cfladd_temp(2))

         !*** vertical momentum diffusion

         imax = maxloc(cfl_vdiffuk_block(:,k))
         cfl_temp = cfl_vdiffuk_block(imax(1),k)
         cfl_vdiffuk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_vdiffuk(k)) then
            cfladd_temp(1:2) = cfladd_vdiffuk_block(:,imax(1),k)
         endif
         cfladd_vdiffuk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_vdiffuk(2,k) = global_maxval(cfladd_temp(2))

         !*** horizontal tracer diffusion

         imax = maxloc(cfl_hdifftk_block(:,k))
         cfl_temp = cfl_hdifftk_block(imax(1),k)
         cfl_hdifftk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_hdifftk(k)) then
            cfladd_temp(1:2) = cfladd_hdifftk_block(:,imax(1),k)
         endif
         cfladd_hdifftk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_hdifftk(2,k) = global_maxval(cfladd_temp(2))

         !*** horizontal momentum diffusion

         imax = maxloc(cfl_hdiffuk_block(:,k))
         cfl_temp = cfl_hdiffuk_block(imax(1),k)
         cfl_hdiffuk(k) = global_maxval(cfl_temp)

         cfladd_temp = 0
         if (cfl_temp == cfl_hdiffuk(k)) then
            cfladd_temp(1:2) = cfladd_hdiffuk_block(:,imax(1),k)
         endif
         cfladd_hdiffuk(1,k) = global_maxval(cfladd_temp(1))
         cfladd_hdiffuk(2,k) = global_maxval(cfladd_temp(2))

      end do

!-----------------------------------------------------------------------
!
!     compute global maximum cfl numbers from level maxima
!
!-----------------------------------------------------------------------

      do k = 1,km

         !*** advective cfl numbers

         if (cfl_advuk(k) > cfl_advu) then
            cfl_advu = cfl_advuk(k)
            cfladd_advu(1) = cfladd_advuk(1,k)
            cfladd_advu(2) = cfladd_advuk(2,k)
            cfladd_advu(3) = k
         endif

         if (cfl_advvk(k) > cfl_advv) then
            cfl_advv = cfl_advvk(k)
            cfladd_advv(1) = cfladd_advvk(1,k)
            cfladd_advv(2) = cfladd_advvk(2,k)
            cfladd_advv(3) = k
         endif

         if (cfl_advwk(k) > cfl_advw) then
            cfl_advw = cfl_advwk(k)
            cfladd_advw(1) = cfladd_advwk(1,k)
            cfladd_advw(2) = cfladd_advwk(2,k)
            cfladd_advw(3) = k
         endif

         if (cfl_advtk(k) > cfl_advt) then
            cfl_advt = cfl_advtk(k)
            cfladd_advt(1) = cfladd_advtk(1,k)
            cfladd_advt(2) = cfladd_advtk(2,k)
            cfladd_advt(3) = k
         endif

         !*** vertical diffusion cfl numbers

         if (cfl_vdifftk(k) > cfl_vdifft) then
            cfl_vdifft = cfl_vdifftk(k)
            cfladd_vdifft(1) = cfladd_vdifftk(1,k)
            cfladd_vdifft(2) = cfladd_vdifftk(2,k)
            cfladd_vdifft(3) = k
         endif

         if (cfl_vdiffuk(k) > cfl_vdiffu) then
            cfl_vdiffu = cfl_vdiffuk(k)
            cfladd_vdiffu(1) = cfladd_vdiffuk(1,k)
            cfladd_vdiffu(2) = cfladd_vdiffuk(2,k)
            cfladd_vdiffu(3) = k
         endif

         !*** horizontal diffusion cfl numbers

         if (cfl_hdifftk(k) > cfl_hdifft) then
            cfl_hdifft = cfl_hdifftk(k)
            cfladd_hdifft(1) = cfladd_hdifftk(1,k)
            cfladd_hdifft(2) = cfladd_hdifftk(2,k)
            cfladd_hdifft(3) = k
         endif

         if (cfl_hdiffuk(k) > cfl_hdiffu) then
            cfl_hdiffu = cfl_hdiffuk(k)
            cfladd_hdiffu(1) = cfladd_hdiffuk(1,k)
            cfladd_hdiffu(2) = cfladd_hdiffuk(2,k)
            cfladd_hdiffu(3) = k
         endif

      end do

!-----------------------------------------------------------------------
!
!     write results to stdout
!
!-----------------------------------------------------------------------

      if (my_task == master_task) then
         write(stdout,blank_fmt)

         name_temp1 = 'zonal advection          '
         write(stdout,cfl_fmt1) name_temp1, cfl_advu, cfladd_advu(:)
         name_temp1 = 'meridional advection     '
         write(stdout,cfl_fmt1) name_temp1, cfl_advv, cfladd_advv(:)
         name_temp1 = 'vertical advection       '
         write(stdout,cfl_fmt1) name_temp1, cfl_advw, cfladd_advw(:)
         name_temp1 = 'total advection          '
         write(stdout,cfl_fmt1) name_temp1, cfl_advt, cfladd_advt(:)
 
         if (cfl_vdifft /= c0) then
            name_temp1 = 'vertical tracer mixing   '
            write(stdout,cfl_fmt1) name_temp1, cfl_vdifft, &
                                               cfladd_vdifft(:)
         endif
         if (cfl_vdifft /= c0) then
            name_temp1 = 'vertical moment. mixing  '
            write(stdout,cfl_fmt1) name_temp1, cfl_vdiffu, &
                                               cfladd_vdiffu(:)
         endif
 
         name_temp1 = 'horizontal tracer mixing '
         write(stdout,cfl_fmt1) name_temp1, cfl_hdifft, cfladd_hdifft(:)
         name_temp1 = 'horiz. momentum mixing   '
         write(stdout,cfl_fmt1) name_temp1, cfl_hdiffu, cfladd_hdiffu(:)
 
         if (cfl_all_levels) then
            write(stdout,'(a27)') 'CFL numbers at all levels: '
            do k=1,km
               name_temp1 = 'total advection          '
               write(stdout,cfl_fmt1) name_temp1, cfl_advtk(k), &
                                      cfladd_advtk(:,k), k
               name_temp1 = 'horizontal tracer mixing '
               write(stdout,cfl_fmt1) name_temp1, cfl_hdifftk(k), &
                                      cfladd_hdifftk(:,k), k
            end do
         endif
 
!-----------------------------------------------------------------------
!
!        write results to output file
!
!-----------------------------------------------------------------------

         open(diag_unit, file=diag_outfile, status='old', &
              position='append')

         name_temp1 = 'cfl_advu'
         name_temp2 = 'cfl_advv'
         write(diag_unit,cfl_fmt2) tday, cfl_advu, name_temp1
         write(diag_unit,cfl_fmt2) tday, cfl_advv, name_temp2
         name_temp1 = 'cfl_advw'
         name_temp2 = 'cfl_advt'
         write(diag_unit,cfl_fmt2) tday, cfl_advw, name_temp1
         write(diag_unit,cfl_fmt2) tday, cfl_advt, name_temp2

         if (cfl_vdifft /= c0) then
            name_temp1 = 'cfl_vdifft'
            name_temp2 = 'cfl_vdiffu'
            write(diag_unit,cfl_fmt2) tday, cfl_vdifft, name_temp1
            write(diag_unit,cfl_fmt2) tday, cfl_vdiffu, name_temp2
         endif

         name_temp1 = 'cfl_hdifft'
         name_temp2 = 'cfl_hdiffu'
         write(diag_unit,cfl_fmt2) tday, cfl_hdifft, name_temp1
         write(diag_unit,cfl_fmt2) tday, cfl_hdiffu, name_temp2

         if (cfl_all_levels) then
            do k=1,km
               write(name_temp1,"('cfl_advtk ',i3)") k
               write(name_temp2,"('cfl_advuk ',i3)") k
               write(diag_unit,cfl_fmt2) tday,cfl_advtk(k),  name_temp1
               write(diag_unit,cfl_fmt2) tday,cfl_hdifftk(k),name_temp2
            end do
         endif

!-----------------------------------------------------------------------
!
!        finished writing max cfl numbers
!
!-----------------------------------------------------------------------

         close(diag_unit)

      endif ! master task

!-----------------------------------------------------------------------
!
!     all done
!
!-----------------------------------------------------------------------

   endif ! ldiag_cfl

!-----------------------------------------------------------------------
!EOC

 end subroutine cfl_check

!***********************************************************************
!BOP
! !IROUTINE: check_KE
! !INTERFACE:

 function check_KE(ke_thresh)

! !DESCRIPTION:
!  Checks the kinetic energy diagnostic to determine whether it
!  exceeds a threshold value.  Also checks for negative values.
!  Exceeding these ranges results in a true value returned.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), intent(in) :: &
      ke_thresh             ! threshold value for kinetic energy

! !OUTPUT PARAMETERS:

   logical (log_kind) :: &
      check_KE           ! true if k.e. exceeds threshold or < 0

!EOP
!BOC
!-----------------------------------------------------------------------

   check_KE = .false.  ! default value
   if (diag_ke > ke_thresh .or. diag_ke < 0) check_KE = .true.

!-----------------------------------------------------------------------
!EOC

 end function check_KE

!***********************************************************************

 end module diagnostics

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
